library(tidyverse)
library(cowplot)
library(patchwork)
library(raster)
library(spatstat)
library(KrigR)
library(lubridate)
library(sf)
library(sp)
library(viridis)
library(ggthemes)
library(terra)
library(MetBrewer)
library(ncdf4)Weather data
Downscaling and extracting weather variables
Packages
White mold data
If you have not downloaded the data necessary to run these analysis, please refer to getting started section before runnig the code below
wm_data = read.csv("data_white-mold/WhiteMoldSurveyWrangledData.csv")Removing missing coordinates
wm_data2 = wm_data %>%
filter(!is.na(latitude)) Names of the counties in the dataset
counties = unique(wm_data2$county)
counties[1] "Genesee" "Niagara" "Orleans" "Livingston" "Wyoming"
[6] "Ontario" "Yates" "Chautauqua" "Monroe"
central_ny = c("wyoming", "livingston","ontario","yates")Ploting the location of each field
#covert the names to lowercase
map_snap_fun = function(yearr){
counties_lc = tolower(counties)
ny_map_data = map_data('county', region = 'new york')
snap_map = ny_map_data %>%
filter(subregion %in% counties_lc) %>%
mutate(region2 = case_when(subregion %in% central_ny ~ "Central lakes",
!subregion %in% central_ny ~ "Great lakes")) %>%
ggplot()+
geom_polygon(data = ny_map_data,
aes(x=long, y = lat, group = group),
fill= "white",
color = "black",
size =0.3)+
geom_polygon(aes(x=long, y = lat, group = group,
# fill=region2
),
fill= "#fafced",
color = "black",
size =0.3)+
geom_point(data = wm_data2 %>%
group_by(subject) %>%
filter(dap == max(dap),
!is.na(wm)) %>%
ungroup() %>%
filter(year == yearr) %>%
arrange(wm),
size = 1,
aes(longitude,latitude, color = wm>0))+
coord_map(xlim = c(-80,-76.7),
ylim = c(41.8, 43.5))+
# scale_fill_manual(values = c("gray85", "gray70"))+
# scale_color_colorblind(labels = c("Absent", "Present"))+
scale_color_manual(labels = c("Absent", "Present"),
values = c("#009E73", "#E69F00"))+
guides(color = guide_legend(override.aes = list(size=3)),
fill = "none")+
theme_minimal()+
labs(x = "Longitude",
y = "Latitude",
fill = "Climate division",
color = "White mold",
title = paste(yearr))+
# facet_wrap(~year, ncol = 1)+
theme(legend.position = "right",
legend.text = element_text(size=7),
axis.title = element_text(size=7),
axis.text = element_text(size=7),
plot.title = element_text(size=10))
snap_map
}ny_map_data = map_data('county', region = 'new york')
counties_lc = tolower(counties)
ny_map = ny_map_data %>%
filter(subregion %in% counties_lc) %>%
mutate(region2 = case_when(subregion %in% central_ny ~ "Central lakes",
!subregion %in% central_ny ~ "Great lakes")) %>%
ggplot()+
geom_polygon(data = ny_map_data,
aes(x=long, y = lat, group = group),
fill= "white",
color = "black",
size =0.3
)+
geom_polygon(aes(x=long, y = lat, group = group,
# fill=region2
),
fill= "#fafced",
color = "black",
size =0.3)+
annotate("rect", xmin = -80, xmax = -76.5,
ymin = 41.8, ymax = 43.5,
color = "black",
size = 0.3,
alpha = 0)+
# scale_size_manual(values = c(0.1,0.3))+
# coord_map(xlim = c(-80,-73),
# ylim = c(40, 45))+
coord_map()+
# scale_fill_manual(values = c("gray85", "gray70"))+
theme_void()+
labs(x = "Longitude",
y = "Latitude",
fill = "Climate division")+
theme(legend.position = c(0.1,0.8))Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
ny_map (map_snap_fun(yearr = "2006") +
map_snap_fun(yearr = "2007")+
map_snap_fun(yearr = "2008")+plot_layout(ncol = 2, guides = 'collect')&
theme(legend.position = "bottom",
legend.justification = c("left","center"))) +
(ny_map+plot_layout(guides = 'keep')) +
# plot_layout(ncol = 2)+
plot_annotation(tag_levels = "A")&
theme(legend.text = element_text(size=6),
legend.title = element_text(size=8),
legend.key.size = unit(0.4, 'cm'),
plot.title = element_text(face = "bold"),
plot.tag = element_text(face = "bold"))ggsave("figs/maps/maps_fields.png", dpi= 900, height = 6, width = 7, bg = "white")
ggsave("figs/maps/maps_fields.pdf", dpi= 900, height = 6, width = 7, bg = "white")ERA5
ERA5-Land hourly data from 1950 to present Description
0.1° x 0.1°; Native resolution is 9 km
Global
Hourly
January 1950 to present
ERA5 hourly data on pressure levels from 1979 to present Description
Reanalysis: 0.25° x 0.25°
Global
Hourly
1979 to present
Agrometeorological indicators from 1979 to present derived from reanalysis Description
0.1° x 0.1°
Global
Daily
1979 to present
There is a package for downloading the data using R: link
ERA5-Land hourly data from 1950 to present
Besides there is a package for downloading data from Era5, I chose using their website to download the data. The process, using their website or the R pakage require making a request and waiting a long time to get the data.
These are the codes for each variable inside the NETCDF files.
u10,d2m,t2m,lai_hv,lai_lv,src,skt,stl2,sp,tp,swvl1,swvl2,swvl3,swvl4
For now, we are going to use only these variables:
t2m: 2m temperatured2m: 2m dew point temperaturesp: Surface pressureswvl1: Soil moisturestl1: Soil temperature
Relative humidity can be calculated from t2m, d2m, and sp.
Importing data
Here I load the data for each data variable and gather into a single raster object.
2m temperature
t2m_all = raster::stack("data_era5/era5_NY_2006-2008.nc", varname = "t2m")2m dewpoint temperature
d2m_all = raster::stack("data_era5/era5_NY_2006-2008.nc", varname = "d2m")Surface pressure
sp_all = raster::stack("data_era5/era5_NY_2006-2008.nc", varname = "sp")soil moisture
swvl1_all = raster::stack("data_era5/era5_NY_2006-2008.nc", varname = "swvl1")soil temperature
stl1_all = raster::stack("data_era5/era5_NY_2006-2008_2.nc", varname = "stl1")Kringing Test
Weather data
Here I load the temperature data, select one layer (time and day) and plot it together with the New York map
r <- raster::stack("data_era5/era5_NY_2006-2008.nc", varname = "t2m")
# a = r$X2008.05.01.00.00.00
a = subset(r,1)
raster::plot(a)
points(wm_data2$longitude, wm_data2$latitude)
maps::map('county', region = 'new york', col = "#5E610B", add = TRUE)length(names(r))[1] 15408
NY shapefile
We will need the New York shapefile to download the digital elevation map, which should be used for kriging. I also filter only the counties that have been surveyed for white mold.
Dir.StateShp <- file.path("data_era5/test")# file.path("shape_files/nc_files")
# ny_shape1 = readOGR("shape_files/cugir-007865/cugir-007865/cty036.shp")
ny_shape1 = sf::read_sf("shape_files/cugir-007865/cugir-007865/cty036.shp")
ny_shape = ny_shape1[ny_shape1$NAME %in% c(unique(wm_data2$county)
# "Yates", "Steuben","Allegany", "Cattaraugus","Erie"
),]
ggplot()+
geom_sf(data =ny_shape )+
geom_point(aes(wm_data2$longitude, wm_data2$latitude))+
theme_nothing()Masking
In this section, we create a mask to filter the dataset, retaining only data for the selected counties. The process involves the following steps:
- Convert the selected county shapes into a simple feature object for visualization.
ny_shape_st = st_as_sfc(ny_shape[1])
plot(ny_shape_st)- Merge the geometries to create a single boundary.
ny_shape_st_nobound = st_union(ny_shape_st)
plot(ny_shape_st_nobound)- Apply a buffer around the selected area to ensure a smooth mask transition.
ny_shape_buffer = st_sf(st_buffer(ny_shape_st_nobound,dist = 10000, max_cells = 500))
plot(ny_shape_buffer)- Use the buffer to mask the raster, setting all selected cells to 1.
a11 = terra::mask(a, ny_shape_buffer, updatevalue=NA, touches =T)
a11_values = values(a11)
a11_values[!is.na(values(a11))] = 1
values(a11) = a11_values
plot(a11)- Convert the masked raster into a polygon. This will serve as the mask for filtering all datasets consistently.
a21 = raster::rasterToPolygons(a11, dissolve=T)- Apply the mask to filter the raster dataset, ensuring only the relevant regions remain.
a2 = crop(mask(a,a21, touches =T),a21)- Visualize the results:
- The first plot shows the initial masked raster.
- The second plot shows the mask polygon.
- The third plot displays the filtered raster dataset.
par(mfrow=c(1,3))
plot(a11, main = "NULL raster file")
plot(a21, main = "Mask")
plot(a2, main = "Masked raster data")ori_data_df = as.data.frame(a2, xy = T)
colnames(ori_data_df)[3] <- "values"
ori_g= ori_data_df %>%
filter(!is.na(values)) %>%
ggplot()+
geom_tile(aes(x, y, fill = values))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
color = "black")+
scale_fill_viridis(option="B")+
coord_sf(xlim = c(-80,-76.8), ylim = c(41.8,43.5))+
geom_point(data = wm_data2,
shape = 21,
color = "white",
fill = "black",
aes(longitude,latitude))+
labs(x = "Longitude",
y = "Latitude",
title= "Era5 original")+
theme_minimal()
ori_gDownloading Digital elevation model (DEM)
Function to plot the DEM
source("functions/Plot_Covs.R")Downloading DEM
# Covs_ls <- download_DEM(
# Train_ras = a2,
# Target_res = .02,
# Shape = ny_shape,
# Dir = Dir.StateShp,
# Keep_Temporary = TRUE
# )
Covs_ls <- CovariateSetup(
Training = as(a2, "SpatRaster"),
Target = .02,
Extent = a21,
Dir = "data_era5",
Keep_Global = TRUE
)###### Downloading global covariate data
[1] "GMTED2010 covariate data already downloaded."
###### Resampling Data
# nc <- nc_open("data_era5/GMTED2010.nc")
# names(nc[['var']])
# Covs_ls = raster::stack(x = "data_era5/GMTED2010.nc", varname = "GMTED2010")
# names(Covs_ls)
# Plot_Covs(Covs_ls)
Plot.Covariates(Covs_ls)Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.
Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.
# class(a2)Visualization
…Using ggplot2
as.data.frame(Covs_ls[[2]], xy = T) %>%
mutate(DEM = GMTED2010) |>
filter(!is.na(DEM)) %>%
ggplot() +
geom_tile(aes(x,y,fill = DEM))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
color = "black")+
geom_point(data = wm_data2, size = 0.1,
aes(longitude,latitude))+
scale_fill_gradientn(colors=met.brewer("Homer2", direction = -1))+
coord_sf(xlim = c(-80,-76.8), ylim = c(41.8,43.5))+
labs(x = "Longitude",
y = "Latitude",
title= "Era5 original")+
theme_void()+
theme(legend.position = "none")ggsave("figs/maps/elevation.png", dpi = 600, height = 5, width = 7, bg = "white")Kringing
This function performs kriging of the weather variables as function of the DEM data
a2_SpatRaster = as(a2, "SpatRaster")
KrigStart <- Sys.time()
State_Krig <- Kriging(
Data =a2_SpatRaster, # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 1, # we want to krig using three cores to speed this process up
FileName = "State_Shape", # save the finished file as this _t2m_2008.nc
Dir = "data_era5/kriging" # where to save the output
)###### Checking your Kriging Specification
[1] "A file with the name State_Shape_Kriged.nc already exists in data_era5/kriging."
[1] "Loading this file for you from the disk."
[1] "A file with the name State_Shape_StDev.nc already exists in data_era5/kriging."
[1] "Loading this file for you from the disk."
KrigStop <- Sys.time()
KrigStop-KrigStartTime difference of 0.116817 secs
Plot.Kriged(State_Krig)# plot(Covs_ls[[1]])Visualization
Krigs = State_Krig$Prediction
# Krigs = State_Krig$Kriging_SE
Krig_df <- as.data.frame(Krigs[[1]], xy = TRUE)
colnames(Krig_df)[3] <- "values"
krig_g = Krig_df %>%
filter(!is.na(values)) %>%
ggplot()+
geom_tile(aes(x, y, fill = values))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
color = "black")+
scale_fill_viridis(option="B")+
# scale_color_viridis(option="B")+
geom_point(data = wm_data2,
shape = 21,
color = "white",
fill = "black",
aes(longitude,latitude))+
coord_sf(xlim = c(-80,-76.8), ylim = c(41.8,43.5))+
labs(x = "Longitude",
y = "Latitude",
title= "Era5 kriged")
ori_g + krig_g +
plot_layout(ncol = 1) &
theme_minimal()Daily summaries
par(mfrow=c(1,3))
plot(mean(t2m_all[[1:24]]))
maps::map('county', region = 'new york', col = "#5E610B", add = TRUE)
plot(min(t2m_all[[1:24]]))
maps::map('county', region = 'new york', col = "#5E610B", add = TRUE)
plot(max(t2m_all[[1:24]]))
maps::map('county', region = 'new york', col = "#5E610B", add = TRUE)# days_in_month(month)
year = 2006:2008
month = 4:10
day = 01:31
hour = as.numeric(seq(0,23,1))
data_era5 = expand_grid(year, month, day) %>%
dplyr::filter(month != 4 | day != 31,
month != 6 | day != 31,
month != 9 | day != 31)%>%
unite(date, year, month, day,sep="-",remove = F ) %>%
mutate(date = as.Date(date))Function for calculating the summaries
mean_raster = function(i, stack_obj){
mean(stack_obj[[i:(i+23)]])
}
min_raster = function(i, stack_obj){
min(stack_obj[[i:(i+23)]])
}
max_raster = function(i, stack_obj){
max(stack_obj[[i:(i+23)]])
}Lapply
days_i = seq(1,length(data_era5$day)*24,by = 24)
time_start = Sys.time()
# dew point
aa1 = lapply(days_i, mean_raster, stack_obj= d2m_all)
d2m_mean_daily_stack = crop(mask(stack(aa1), a21),a21)
writeCDF(rast(d2m_mean_daily_stack),
"data_era5/daily_summaries/d2m_mean_daily.nc",
overwrite=TRUE)
# temperature
### mean
aa2 = lapply(days_i, mean_raster, stack_obj= t2m_all)
t2m_mean_daily_stack = crop(mask(stack(aa2), a21),a21)
writeCDF(rast(t2m_mean_daily_stack),
"data_era5/daily_summaries/t2m_mean_daily.nc",
overwrite=TRUE)
### minimum
aa3 = lapply(days_i, min_raster, stack_obj= t2m_all)
t2m_min_daily_stack = crop(mask(stack(aa3), a21),a21)
writeCDF(rast(t2m_min_daily_stack),
"data_era5/daily_summaries/t2m_min_daily.nc",
overwrite=TRUE)
### maximum
aa4 = lapply(days_i, max_raster, stack_obj= t2m_all)
t2m_max_daily_stack = crop(mask(stack(aa4), a21),a21)
writeCDF(rast(t2m_max_daily_stack),
"data_era5/daily_summaries/t2m_max_daily.nc",
overwrite=TRUE)
# Surface pressure
aa5 = lapply(days_i, mean_raster, stack_obj= sp_all)
sp_mean_daily_stack = crop(mask(stack(aa5), a21),a21)
writeCDF(rast(sp_mean_daily_stack),
"data_era5/daily_summaries/sp_mean_daily.nc",
overwrite=TRUE)
# Soil moisture
aa6 = lapply(days_i, mean_raster, stack_obj= swvl1_all)
swvl1_mean_daily_stack = crop(mask(stack(aa6), a21),a21)
writeCDF(rast(swvl1_mean_daily_stack),
"data_era5/daily_summaries/swvl1_mean_daily.nc",
overwrite=TRUE)
# soil temperature
aa7 = lapply(days_i, mean_raster, stack_obj= stl1_all)
stl1_mean_daily_stack = crop(mask(stack(aa7), a21),a21)
writeCDF(rast(stl1_mean_daily_stack),
"data_era5/daily_summaries/stl1_mean_daily.nc",
overwrite=TRUE)
time_end = Sys.time()
time_end -time_startplot(d2m_mean_daily_stack$X1)
# plot(stl1_mean_daily_stack$X1)
maps::map('county', region = 'new york', add = TRUE)Download DEM
Kriging the daily summaries
Dir.krigd_var <- file.path("data_era5/kriged")Dew point
d2m_mean_daily_Krig = Kriging(
Data = as(d2m_mean_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "d2m_mean_daily_krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)Mean Temeperature
t2m_mean_daily_Krig = Kriging(
Data = as(t2m_mean_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "t2m_mean_daily_Krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)Minimum temperature
t2m_min_daily_Krig = Kriging(
Data = as(t2m_min_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "t2m_min_daily_Krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)Maximum temperature
t2m_max_daily_Krig = Kriging(
Data = as(t2m_max_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "t2m_max_daily_Krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)Surface pressure
sp_mean_daily_Krig = Kriging(
Data = as(sp_mean_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "sp_mean_daily_Krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)Soil temperature
stl1_mean_daily_Krig = Kriging(
Data = as(stl1_mean_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "stl1_mean_daily_Krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)Soil moisture
swvl1_mean_daily_Krig = Kriging(
Data = as(swvl1_mean_daily_stack, "SpatRaster"), # what to krig
Covariates_training = Covs_ls[[1]], # covariates at training resolution
Covariates_target = Covs_ls[[2]], # covariates at target resolution
Equation = "GMTED2010", # the covariate(s) we want to use
Keep_Temporary = FALSE, # delete temporary krigs of layers
nmax = 40, # degree of localisation
Cores = 5, # we want to krig using three cores to speed this process up
FileName = "swvl1_mean_daily_Krig.nc", # save the finished file as this _t2m_2008.nc
Dir = Dir.krigd_var#"data_era5/kriging" # where to save the output
)d2m_mean_daily_Krig = stack("data_era5/kriged/d2m_mean_daily_Krig.nc")
t2m_mean_daily_Krig = stack("data_era5/kriged/t2m_mean_daily_Krig.nc")
t2m_min_daily_Krig = stack("data_era5/kriged/t2m_min_daily_Krig.nc")
t2m_max_daily_Krig = stack("data_era5/kriged/t2m_max_daily_Krig.nc")
sp_mean_daily_Krig = stack("data_era5/kriged/sp_mean_daily_Krig.nc")
swvl1_mean_daily_Krig = stack("data_era5/kriged/swvl1_mean_daily_Krig.nc")
stl1_mean_daily_Krig = stack("data_era5/kriged/stl1_mean_daily_Krig.nc")Ploting the kingued maps (first layer)
par(mfrow=c(2,4))
plot(t2m_mean_daily_Krig$X1, main = "Mean Temp")
plot(t2m_min_daily_Krig$X1, main = "Min Temp")
plot(t2m_max_daily_Krig$X1, main = "Max Temp")
plot(d2m_mean_daily_Krig$X1, main = "Dew Point")
plot(sp_mean_daily_Krig$X1, main = "Surface Pressure")
plot(swvl1_mean_daily_Krig$X1, main = "Soil Moisture")
plot(stl1_mean_daily_Krig$X1, main = "Soil Temp")Using ggplot
# t2m_max_daily_stack
Krigs = t2m_max_daily_stack$X1
Krig_df <- as.data.frame(Krigs[[1]], xy = TRUE)
colnames(Krig_df)[3] <- "values"
max_t_ori = Krig_df %>%
filter(!is.na(values)) %>%
ggplot()+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= "white",
color = "black")+
geom_tile(aes(x, y, fill = values-273.15))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
color = "black")+
scale_fill_viridis(option="B", limits = c(5,18),breaks =seq(5, 18, by = 3))+
# scale_color_viridis(option="B")+
geom_point(data = wm_data2,
shape = 21,
color = "white",
fill = "black",
aes(longitude,latitude))+
coord_sf(xlim = c(-80,-76.8), ylim = c(41.8,43.5))+
labs(x = "Longitude",
y = "Latitude",
fill = "Max Temperature (°C)",
title= "Era5 Native Resolution (0.1° x 0.1°)")
# max_t_oriKrigs = t2m_max_daily_Krig$X1
Krig_df <- as.data.frame(Krigs[[1]], xy = TRUE)
colnames(Krig_df)[3] <- "values"
max_t_krig = Krig_df %>%
filter(!is.na(values)) %>%
ggplot()+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= "white",
color = "black")+
geom_tile(aes(x, y, fill = values-273.15))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
color = "black")+
scale_fill_viridis(option="B", limits = c(5,18), breaks =seq(5, 18, by = 3))+
# scale_color_viridis(option="B")+
geom_point(data = wm_data2,
shape = 21,
color = "white",
fill = "black",
aes(longitude,latitude))+
coord_sf(xlim = c(-80,-76.8), ylim = c(41.8,43.5))+
labs(x = "Longitude",
y = "Latitude",
fill = "Max Temperature (°C)",
title= "Era5 Kriged (0.02° x 0.02°)")
max_t_ori + max_t_krig +
plot_layout(ncol = 2, guides = "collect") &
theme_map()&
theme(legend.position = "bottom",
legend.text = element_text(size = 9))# ggsave("figs/max_t.png",dpi=900, height = 8, width=14, bg = "white")Relative humidity
How to calculate relative humidity here link and here link
\[RH = \frac{es(d2m)}{es(t2m)}*100\]
\(es()\) is Saturation vapor pressure and can be calculated as
\[es(T) = \alpha_1*exp( \alpha_3*(\frac{T-t0}{T-\alpha_4}))\],
where \(T\) is the temperature, \(\alpha_1 = 611.21\), \(\alpha_2 = 17.502\), \(\alpha_3 = 32.19\).
Using these resources, I build the function to calculates RH from Temperature and dew point.
es =function(Temp){
t0 = 273.16
alpha1 = 611.21
alpha3 = 17.502
alpha4 = 32.19
alpha1*exp( alpha3*((Temp-t0)/(Temp-alpha4)))
}
#test
Tdp = 14+273.15
T2m = 31+273.15
es(Tdp)/es(T2m)*100[1] 35.55801
Creating the raster object of relative humidity
layer_i = 1:nlayers(d2m_mean_daily_Krig)
calculate_RH = function(i, d2m_stack, t2m_stack){
es(d2m_stack[[i]])/es(t2m_stack[[i]])*100
}
list_rh = lapply(layer_i, calculate_RH, d2m_stack = d2m_mean_daily_Krig,t2m_stack= t2m_mean_daily_Krig)
rh_mean_daily_Krig = brick(list_rh)
writeCDF(rast(rh_mean_daily_Krig),
"data_era5/kriged/rh_mean_daily_Krig.nc",
overwrite=TRUE)rh_mean_daily_Krig = stack("data_era5/kriged/rh_mean_daily_Krig.nc")
plot(rh_mean_daily_Krig$X1, main = "Relative Humidity")Extracting data
wm_data2_uni = wm_data2 %>%
group_by(subject) %>%
slice(1L)coords<-data.frame(lon=wm_data2_uni$longitude, lat=wm_data2_uni$latitude)
coordinates(coords)<-c("lon","lat")- Dew point
ext_d2m = extract(d2m_mean_daily_Krig, coords)
colnames(ext_d2m) = as.character(data_era5$date)
ext_d2m_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_d2m )
d2m_mean_wm = as.data.frame(ext_d2m_coord) %>%
mutate(subject =wm_data2_uni$subject) %>%
pivot_longer(3:644,
values_to = "d2m",
names_to = "date")%>%
mutate(date = as.Date(date))
d2m_mean_wm# A tibble: 241,392 × 5
longitude latitude subject date d2m
<dbl> <dbl> <int> <date> <dbl>
1 -77.9 42.9 1 2006-04-01 282.
2 -77.9 42.9 1 2006-04-02 271.
3 -77.9 42.9 1 2006-04-03 276.
4 -77.9 42.9 1 2006-04-04 270.
5 -77.9 42.9 1 2006-04-05 266.
6 -77.9 42.9 1 2006-04-06 272.
7 -77.9 42.9 1 2006-04-07 275.
8 -77.9 42.9 1 2006-04-08 268.
9 -77.9 42.9 1 2006-04-09 268.
10 -77.9 42.9 1 2006-04-10 270.
# ℹ 241,382 more rows
- Mean Temperature
ext_t2m_mean = extract(t2m_mean_daily_Krig, coords)
colnames(ext_t2m_mean) = as.character(data_era5$date)
ext_t2m_mean_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_t2m_mean)
t2m_mean_wm = as.data.frame(ext_t2m_mean_coord) %>%
mutate(subject =wm_data2_uni$subject)%>%
pivot_longer(3:644,
values_to = "t2m_mean",
names_to = "date")%>%
mutate(date = as.Date(date))- Max Temperature
ext_t2m_max = extract(t2m_max_daily_Krig, coords)
colnames(ext_t2m_max) = as.character(data_era5$date)
ext_t2m_max_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_t2m_max)
t2m_max_wm = as.data.frame(ext_t2m_max_coord) %>%
mutate(subject =wm_data2_uni$subject)%>%
pivot_longer(3:644,
values_to = "t2m_max",
names_to = "date")%>%
mutate(date = as.Date(date))- Min Temperature
ext_t2m_min = extract(t2m_min_daily_Krig, coords)
colnames(ext_t2m_min) = as.character(data_era5$date)
ext_t2m_min_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_t2m_min)
t2m_min_wm = as.data.frame(ext_t2m_min_coord) %>%
mutate(subject =wm_data2_uni$subject)%>%
pivot_longer(3:644,
values_to = "t2m_min",
names_to = "date")%>%
mutate(date = as.Date(date))- Pressure
ext_sp = extract(sp_mean_daily_Krig, coords)
colnames(ext_sp) = as.character(data_era5$date)
ext_sp_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_sp )
sp_mean_wm = as.data.frame(ext_sp_coord) %>%
mutate(subject =wm_data2_uni$subject)%>%
pivot_longer(3:644,
values_to = "sp",
names_to = "date")%>%
mutate(date = as.Date(date))- Soil moisture
ext_sm = extract(swvl1_mean_daily_Krig, coords)
colnames(ext_sm) = as.character(data_era5$date)
ext_sm_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_sm )
sm_mean_wm = as.data.frame(ext_sm_coord)%>%
mutate(subject =wm_data2_uni$subject) %>%
pivot_longer(3:644,
values_to = "sm",
names_to = "date")%>%
mutate(date = as.Date(date))- Relative humidity
ext_rh = extract(rh_mean_daily_Krig, coords)
colnames(ext_rh) = as.character(data_era5$date)
ext_rh_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_rh )
rh_mean_wm = as.data.frame(ext_rh_coord)%>%
mutate(subject =wm_data2_uni$subject)%>%
pivot_longer(3:644,
values_to = "rh",
names_to = "date") %>%
mutate(date = as.Date(date))- Soil temperature
ext_st = extract(stl1_mean_daily_Krig, coords)
colnames(ext_st) = as.character(data_era5$date)
ext_st_coord = cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude,
ext_st )
st_mean_wm = as.data.frame(ext_st_coord) %>%
mutate(subject =wm_data2_uni$subject)%>%
pivot_longer(3:644,
values_to = "st",
names_to = "date") %>%
mutate(date = as.Date(date))Binding data sets
weather_all_wm = d2m_mean_wm %>%
bind_cols(t2m_mean_wm[,5],
t2m_max_wm[,5],
t2m_min_wm[,5],
sp_mean_wm[,5],
sm_mean_wm[,5],
rh_mean_wm[,5],
st_mean_wm[,5])
weather_all_wm# A tibble: 241,392 × 12
longitude latitude subject date d2m t2m_mean t2m_max t2m_min sp
<dbl> <dbl> <int> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
1 -77.9 42.9 1 2006-04-01 282. 284. 289. 281. 97567.
2 -77.9 42.9 1 2006-04-02 271. 279. 285. 273. 98578.
3 -77.9 42.9 1 2006-04-03 276. 282. 290. 277. 97387.
4 -77.9 42.9 1 2006-04-04 270. 277. 284. 273. 97250.
5 -77.9 42.9 1 2006-04-05 266. 273. 277. 270. 97416.
6 -77.9 42.9 1 2006-04-06 272. 276. 283. 273. 97868.
7 -77.9 42.9 1 2006-04-07 275. 278. 281. 275. 97193.
8 -77.9 42.9 1 2006-04-08 268. 274. 277. 271. 97920.
9 -77.9 42.9 1 2006-04-09 268. 275. 282. 270. 98859.
10 -77.9 42.9 1 2006-04-10 270. 278. 288. 272. 99093.
# ℹ 241,382 more rows
# ℹ 3 more variables: sm <dbl>, rh <dbl>, st <dbl>
Saving data
data_full = weather_all_wm %>%
dplyr::select(-latitude, -longitude) %>%
full_join(wm_data2, by = "subject")
write.csv(data_full, "data_white-mold/data_model_plus_weather.csv",row.names = F)data_full = read.csv("data_white-mold/data_model_plus_weather.csv") %>%
mutate(date = as.Date(date),
sampling.date = as.Date(sampling.date),
planting.date = as.Date(planting.date))filter between planting and sampling dates
data_fill_filtered = data_full %>%
filter(date >= (planting.date-30) & date <= sampling.date)write.csv(data_fill_filtered, "data_white-mold/data_model_plus_weather_filtered.csv",row.names = F)data_fill_filtered = read.csv("data_white-mold/data_model_plus_weather_filtered.csv")Session Info
sessionInfo()R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=Portuguese_Brazil.utf8 LC_CTYPE=Portuguese_Brazil.utf8
[3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C
[5] LC_TIME=Portuguese_Brazil.utf8
time zone: America/Sao_Paulo
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] ncdf4_1.23 MetBrewer_0.2.0 terra_1.8-5
[4] ggthemes_5.1.0 viridis_0.6.5 viridisLite_0.4.2
[7] sf_1.0-19 KrigR_0.9.4 spatstat_3.3-0
[10] spatstat.linnet_3.2-3 spatstat.model_3.3-3 rpart_4.1.23
[13] spatstat.explore_3.3-3 nlme_3.1-164 spatstat.random_3.3-2
[16] spatstat.geom_3.3-4 spatstat.univar_3.1-1 spatstat.data_3.1-4
[19] raster_3.6-30 sp_2.1-4 patchwork_1.3.0
[22] cowplot_1.1.3 lubridate_1.9.3 forcats_1.0.0
[25] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
[28] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[31] ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] DBI_1.2.3 pbapply_1.7-2 deldir_2.0-4
[4] gridExtra_2.3 s2_1.1.7 rlang_1.1.4
[7] magrittr_2.0.3 e1071_1.7-16 compiler_4.4.1
[10] mgcv_1.9-1 systemfonts_1.1.0 maps_3.4.2.1
[13] vctrs_0.6.5 ecmwfr_2.0.3 wk_0.9.4
[16] crayon_1.5.3 pkgconfig_2.0.3 fastmap_1.2.0
[19] labeling_0.4.3 utf8_1.2.4 rmarkdown_2.28
[22] tzdb_0.4.0 ragg_1.3.3 automap_1.1-12
[25] xfun_0.48 cachem_1.1.0 jsonlite_1.8.9
[28] progress_1.2.3 goftest_1.2-3 reshape_0.8.9
[31] spatstat.utils_3.1-1 prettyunits_1.2.0 parallel_4.4.1
[34] R6_2.5.1 stringi_1.8.4 stars_0.6-7
[37] Rcpp_1.0.13 iterators_1.0.14 knitr_1.48
[40] tensor_1.5 zoo_1.8-12 snow_0.4-4
[43] FNN_1.1.4.1 Matrix_1.7-0 splines_4.4.1
[46] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.17.0
[49] abind_1.4-8 yaml_2.3.10 codetools_0.2-20
[52] lattice_0.22-6 intervals_0.15.5 plyr_1.8.9
[55] withr_3.0.2 evaluate_1.0.1 units_0.8-5
[58] proxy_0.4-27 polyclip_1.10-7 xts_0.14.1
[61] pillar_1.9.0 KernSmooth_2.23-24 renv_1.1.2
[64] foreach_1.5.2 generics_0.1.3 spacetime_1.3-2
[67] hms_1.1.3 munsell_0.5.1 scales_1.3.0
[70] class_7.3-22 glue_1.8.0 mapproj_1.2.11
[73] tools_4.4.1 grid_4.4.1 colorspace_2.1-1
[76] cli_3.6.3 spatstat.sparse_3.1-0 gstat_2.1-2
[79] textshaping_0.4.0 fansi_1.0.6 doSNOW_1.0.20
[82] gtable_0.3.5 digest_0.6.37 classInt_0.4-10
[85] htmlwidgets_1.6.4 farver_2.1.2 memoise_2.0.1
[88] htmltools_0.5.8.1 lifecycle_1.0.4 httr_1.4.7